In this notebook we will implement the calculation of CONSTANT HEIGHT STM IMAGES.

The equations

Using the Tersoff-Hammann aproximation (https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.50.1998), the tunneling current is proportional to the local density of states at a certain point close to the material:

$$ I(\vec{r}, V) \propto \int_0^V LDOS(\vec{r}, E_f - V + \epsilon) \,d\epsilon\ $$

The LDOS at a certain energy is simply described as the electron density, computed as $\psi(\vec{r}) \psi^*(\vec{r})$, considering only wavefunctions at that energy:

$$ LDOS(r, E) = \sum_{i} n(E, \epsilon_i) * \psi_i(\vec{r}) \psi_i^*(\vec{r}) $$

Where $n(E, \epsilon)$ is the occupation function that determines the contribution of a given wavefunction. $E$ is the energy for which we want the LDOS and $\epsilon_i$ is the energy of the wavefunction. One possible occupation function is Dirac's delta:

$$ n(E, \epsilon) = \delta(E - \epsilon) = \begin{cases} 1, & \text{if}\ E = \epsilon \\ 0, & \text{otherwise}\ \end{cases}\,. $$

In summary, to compute constant height STM images. We need to:

• Get a Hamiltonian (we will use SIESTA for this).
• Compute the eigenstates (wavefunction coefficients) and eigenvalues (energies).
• Implement a function that computes the STM current $I(\vec{r}, V)$ for any point in space $\vec{r}$ at a given voltage $V$. The full equation is:

$$ I(\vec{r}, V) \propto \int_0^V \sum_{i} \left[n(-V + E, \epsilon_i) * \psi_i(\vec{r}) \psi_i^*(\vec{r})\right] \,dE\ $$

• Take the points of space of a given plane (at constant height), as an STM tip would do.

In this first case, we will compute STM images of a benzene molecule.

Tips

To build the molecule:

You can use the ase (Atomic Simulation Environment) python package. It has a function ase.build.molecule that can generate a benzene molecule. Look at its documentation (https://wiki.fysik.dtu.dk/ase/ase/build/build.html#ase.build.molecule) to understand how to use it. Add at least 5 Ang of vacuum!

To write the molecule to a file:

ase will give you an ase.Atoms object. You can convert it to sisl with sisl.Geometry.new(ase_atoms). Then you can write it to an fdf file as usual.

Plot the molecule before writing it to make sure that it is what you expect.

To run the SIESTA calculation:

You will need the Hamiltonian, so don't forget to tell siesta to store it in a .TSHS file by using TS.HS.Save true.

To read the Hamiltonian:

Once the SIESTA run has finished, you can do:

sisl.get_sile("path/to/my/file.fdf").read_hamiltonian()

as we have done in previous lab sessions.

To compute wavefunctions on the grid:

To compute $\psi (\vec{r})$ you need three things.

1. The eigenstate coefficients. Once you have a hamiltonian, you can get all of them with H.eigenstate(). Then you can loop this object to get each individual eigenstate. Each eigenstate has its energy stored under the .eig property.
2. A grid of points in space. You can create one with sisl.Grid(geometry, shape=(100, 100, 100)). This will create a grid of $100x100x100$ points within the cell of your geometry.
3. A function to project the wavefunction into the grid. The eigenstate object has a wavefunction method (docs) that will project the wavefunction into an already initialized grid.

Remember that we have done this already on the first SIESTA lab, with a water molecule.

To plot the STM images at a certain height:

A function is provided to plot molecules at multiple heights. You can use it as it is or modify it as you wish.

Note: If you have time, you can modify this function to create a plots at multiple voltages. Even a grid of plots at multiple voltages and multiple heights!

To check that things make sense:

You can plot all the eigenstates of your system to try to understand if the STM images make sense with something like:

H.plot.pdos(data_Erange=[-10, 10], nE=1200, Erange=[-10, 10], kgrid=[1,1,1]).split_DOS(on="species+l+m")

Helper: Grids have a fill method to set all the values of the grid to a given value, if that is needed.

Let's get to work

From now on, you have to implement the STM images!

Look at the equation for current and try to think how that can be translated into code. If you get stuck, experiment with whatever ideas you have! There's no harm in making python error, you won't break it ;)

Some imports:

import sisl
import sisl.viz
import numpy as np
import plotly.io as pio
pio.renderers.default="notebook"
from pathlib import Path

#En aquesta part del codi tenim la construcció de la molècula de benzè amb un "vacuum" de valor 5 Ang.

from ase.build import molecule
atoms = molecule('C6H6',5)

#Definim la molècua com "atoms" i, seguidament, creem un arxiu "geom.fdf" en una carpeta que contingui l'informació de la geometria de la molècula.

geometry = sisl.Geometry.new(atoms)
geometry.write("27M/geom.fdf")

#Fem el "plot" de la geometria del benzè per comprovar que és la correcta.

geometry.plot(axes="xy")

#A continuació, definim les variables "H" i "geom" que correspondran al Hamiltonià i a la geometria obtinguda a partir de fer la simulació Siesta.

H=sisl.get_sile("27M/RUN.fdf").read_hamiltonian()
geom=sisl.get_sile("27M/RUN.fdf").read_geometry()

#La següent part consisteix en el càlcul de la intensitat STM en funció de V, que correspona a l'integral de la LDOS. 

#Primerament, definim la funció "LDOS_grid" que depèn del Hamiltonià "H" i del voltatge que voldrem aplicar.

def LDOS_grid(H,V):
    grid=sisl.Grid((100, 100, 100), geometry=geom) #Aquest codi ens permet crear un "grid" final on tindrem aplicats els càlculs.
    eigenstate=H.eigenstate() #Aquí definim els estats propis a partir del Hamiltonià. 
    gridbuit=sisl.Grid((100, 100, 100), geometry=geom) #Amb aquest altre codi el que busquem és crear un grid temporal per cada funció de densitat del benzè.

#En aquesta part, per valors positius i negatius, respectivament, de voltatge aplicat, fem un "loop" per obtenir un "grid" final amb la probabilitat de densitat.
    
    if V>0:
        for eig in eigenstate:
            if eig.eig<=V and eig.eig>0:
                gridbuit.fill(0) #Fem que el "grid" temporal es buidi.
                wf=eig.wavefunction(gridbuit) #Definim la funció d'ona dintre del "grid" buit.
                densitat=gridbuit*gridbuit #Multipliquem el "grid" buit per ell mateix obtenir la seva densita.
                grid=grid+densitat #Obtenim el "grid" final al sumar-li la densitat obtinguda anteriorment.
    if V<0:
        for eig in eigenstate:
            if eig.eig>=V and eig.eig<0:
                gridbuit.fill(0)
                wf=eig.wavefunction(gridbuit)
                densitat=gridbuit*gridbuit
                grid=grid+densitat
    return grid

#Ara, he fet el "plot" del resultat per un valor negatiu i positiu de voltatge per comprovar que funciona.

LDOS_grid(H,-15).plot(axes='xy').show()
LDOS_grid(H,15).plot(axes='xy').show()

The helper function to plot STM images once you have the LDOS grid:

# The helper function
def plot_STM_images(
    LDOS_grid: sisl.Grid, 
    min_height: int = 0, 
    max_height: int = 2,
    steps: int = 9,
    crange = None,
    colorscale: str = None
):
    """Plots constant height STM images at multiple heights

    Parameters
    ----------
    LDOS_grid:
        A grid containing the LDOS corresponding to the voltage that you
        want to plot, for ALL SPACE.
    min_height:
        The height of the first image.
    max_height:
        The height of the last image.
    steps:
        The number of steps between the first and the last image.
    log:
        Whether to plot the log of the values.
    crange:
        The range of the colorscale. E.g. [0, 2]. If None, it is computed
        from the minimum and maximum values of the data.
    colorscale:
        The plotly colorscale to use.

    Examples
    ----------

    This function should be used like:
    >>> grid = ...compute the LDOS in the grid.
    >>> plot_STM_images(grid)
    """

    # Determine all the heights that the user wants to plot
    heights = np.linspace(min_height, max_height, steps)
    # Get the position of the benzene molecule
    z_0 = LDOS_grid.geometry.xyz[:, 2].max()
    
    # Initialize a list of plots
    plots = []
    coloraxes = []
    # Loop through the heights, and for each of them we will create a plot.
    for height in heights:
        # Get the z for which we want the LDOS (benzene position + height)
        z = z_0 + height
        
        # Get the plot for this height
        plot = LDOS_grid.plot(axes="xy", z_range=[z - 0.1, z + 0.1], crange=crange, colorscale=colorscale)

        # And append it to the list
        plots.append(plot)

    # Merge all the plots, with some extra arguments to beautify the plot.
    return sisl.viz.merge_plots(
        *plots, 
        composite_method="subplots", 
        arrange="square", 
        subplot_titles=[f"Height = {height} Ang" for height in heights],
        horizontal_spacing=0,
        vertical_spacing=0.05
    ).update_xaxes(
        visible=False
    ).update_yaxes(
        visible=False
    ).update_layout(
        height=900, title=f"Constant height STM images",
    )

H.plot.pdos(data_Erange=[-10, 10], nE=1200, Erange=[-10, 10], kgrid=[1,1,1]).split_DOS(on="species+l+m")

V=[-10,-1,3.2,7.6]
for i in V:
    plot_STM_images(LDOS_grid(H,i)).update_layout(title="Constant height STM images for "+str(i)+" V").show()

Finalment, comentarem la raó de les diferències que observem en la simulació de les imatges de STM per diferents voltatges aplicats, relacionant-ho amb el plot de PDOS.

-10 V: en aquest cas, tindrem la contribució de tots els estats de la molècula partint del 0, obtenint una representació de la molècula real de benzè molt acurada.

-1 V: en aquest voltatge no tenim cap pic de densitats, obtenint com a resultat una representació nul·la. Ens serveix per confirmar que els resultats tenen sentit.

3.2 V: aquí obtenim una representació que correspon a l'observació dels àtoms de carboni, majoritàriament. És un resultat esperat si mirem el gràfoc de PDOS, ja que tenim un pic verd molt alt que correspon a l'espècie C amb l=1 i m=0.

8 V: finalment, en aquest exemple veiem que tenim una gran contribució del orbital s de l'hidrogen, obtenint una representació on veiem de manera clara els seus diferents orbitals.

A més, cal dir que a mesura que augmentem la distància a la que representem, veurem que l'imatge ees va tornant més difusa, tot i que veurem de manera més marcada els orbitals pz del carboni. Això és degut al fet de que aquests orbitals s'orienten de forma perpendicular al pla de la molècula, trobant-se a una distància més propera a la de visualització.

Bonus

En aquest apartat, el que podem fer és definir una variable que sigui la resta de dos grids amb el rang de voltatge que volem ( en aquest cas, de -5.4 fins -3) i, seguidament, representem el resultat.

Bonus=LDOS_grid(H,-5.4)-LDOS_grid(H,-3)

Bonus.plot(axes='xy').show()

En aquest cas, obtenim la representació dels pics que obtenim de la PDOS en el rang establert. Veiem una gran contribució dels orbitals del carboni px, py i pz i, a més, una forta contribució del orbital s de l'hidrogen.